Appendix B Summary Statement: Agricultural Insurance Loss and Principle Components Analysis

Appendix B documents the exploratory data analysis process, examining the relationships of agricultural commodity loss, at a county level, from 1989-2015, across the three state region of the Pacific Northwest (Washington, Idaho, and Oregon), and then focus in on the 24 county region of the Inland Pacific Northwest (IPNW).

Here we perform a variety of Principle Components Analyses (PCA) to better understand the combined effects of differing damage causes, commodities, counties, and years on overall loss.



Step 1. Data Preparation. Here we load our agricultural insurance loss data and prepare it for PCA Analysis.

Step 2. Data Preparation. Original Insurance Loss Dataset - Pacific Northwest

Step 3. Data Preparation. Aggregated Insurance Loss Dataset - Pacific Northwest

Step 4. 24-county inland Pacific Northwest (iPNW) Study Area.

Step 5. PCA: PNW insurance loss, by COUNTY with COMMODITY loadings: 2001 to 2015

Step 6. PCA: PNW insurance loss, by YEAR with DAMAGE CAUSE loadings: 2001 to 2015

Step 7. PCA: PNW insurance loss, by MONTH with DAMAGE CAUSE loadings: 2001 to 2015

Step 8. PCA: PNW WHEAT insurance loss, by YEAR with DAMAGE CAUSE loadings: 2001 to 2015

Step 9. PCA: PNW WHEAT insurance loss, by COUNTY with DAMAGE CAUSE loadings: 2001 to 2015

Step 10. PCA: PNW APPLES insurance loss, by COUNTY with DAMAGE CAUSE loadings: 2001 to 2015

Step 11. PCA: PNW BARLEY insurance loss, by COUNTY with DAMAGE CAUSE loadings: 2001 to 2015

Step 12. PCA: IPNW insurance loss by COUNTY with DAMAGE CAUSE loadings: 2001 to 2015

Step 13. PCA: IPNW insurance loss, by COUNTY with DAMAGE CAUSE loadings: 2001 to 2015

Step 14. PCA: IPNW insurance loss, by YEAR with DAMAGE CAUSE loadings: 2001 to 2015

Step 15. PCA: IPNW insurance loss, by COUNTY with DAMAGE CAUSE loadings: 2015

Step 16. PCA: IPNW WHEAT loss by YEAR with DAMAGE CAUSE loadings: Temporal Mapping of PC1 and PC2

Step 17. PCA: IPNW WHEAT loss by COUNTY with DAMAGE CAUSE loadings: Spatial Mapping of PC1 and PC2

Step 1: Data Preparation

The first step was to load Pacific Northwest agricultural insurance commodity loss data, aggregate by county, damage cause and commodity, remove zeros, transform the loss values ($) into cube root and logrithmic outputs, and scale and center the data.

Original data file for all insurance claims for the Pacific Northwest, can be found here:

https://github.com/erichseamon/paper-PHD1/blob/master/data/RMA_PNW_2001_2015.csv

Aggregated data file, which summarizes data by year, county, commodity, and damage cause, for just the Pacific Northwest, can be found here:

https://github.com/erichseamon/paper-PHD1/blob/master/data/RMA_damage_PNW_2001_2015.csv

All code that generates the following graphs and tables, can be found here:

https://github.com/erichseamon/paper-PHD1/

library(tab)
library(car)
library(RCurl)
library(lme4)
library(ez)
library(lattice)
library(ggplot2)
library(coefplot2)
library(broom)
library(htmlTable)
library(gridExtra)
library(kableExtra)
library(ggfortify)


#options(scipen=999)

#load original data as txt file 

rma1988 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/1988.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma1989 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/1989.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma1990 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/1990.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma1991 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/1991.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma1992 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/1992.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma1993 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/1993.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma1994 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/1994.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma1995 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/1995.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma1996 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/1996.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma1997 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/1997.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma1998 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/1998.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma1999 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/1999.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2000 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2000.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2001 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2001.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2002 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2002.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2003 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2003.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2004 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2004.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2005 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2005.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2006 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2006.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2007 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2007.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2008 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2008.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2009 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2009.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2010 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2010.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2011 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2011.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2012 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2012.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2013 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2013.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2014 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2014.txt", sep = "|", header = FALSE, strip.white = TRUE)
rma2015 <- read.csv("/dmine/data/USDA/crop_indemnity_txt/2015.txt", sep = "|", header = FALSE, strip.white = TRUE)

#load original data as from github as a joint file - 1989 to 2001 and 2001 to 2015

#RMA_2015 <- read.csv("https://raw.githubusercontent.com/erichseamon/paper-PHD1/master/data/RMA_csv/US_RMA_2001_2015.csv", sep = ",", header = FALSE, strip.white = TRUE)

#RMA_1989 <- read.csv("https://raw.githubusercontent.com/erichseamon/paper-PHD1/master/data/RMA_csv/US_RMA_1989_2000.csv", sep = ",", header = FALSE, strip.white = TRUE)


RMA_1989 <- rbind(rma1989, rma1990, rma1991, rma1992, rma1993, rma1994, rma1995, rma1996, rma1997, rma1998, rma1999, rma2000)

RMA_1989 <- data.frame(RMA_1989[,1], RMA_1989[,3], RMA_1989[,5], RMA_1989[,7], RMA_1989[,12], RMA_1989[,13], RMA_1989[,14], RMA_1989[,15])

colnames(RMA_1989) <- c("year", "state", "county", "commodity", "damagecause", "monthcode", "month", "loss")
#write.csv(RMA_1989, file = "/dmine/code/git/paper-PHD1/data/RMA_csv/US_RMA_1989_2000.csv")
#save(RMA_1989,file="/dmine/code/git/paper-PHD1/data/RMA_Rda/US_RMA_1989_2000.Rda")

RMA <- rbind(rma2001, rma2002, rma2003, rma2004, rma2005, rma2006, rma2007, rma2008, rma2009, rma2010, rma2011, rma2012, rma2013, rma2014, rma2015)

RMA <- data.frame(RMA[,1], RMA[,3], RMA[,5], RMA[,7], RMA[,12], RMA[,13], RMA[,14], RMA[,15], RMA[,16])

colnames(RMA) <- c("year", "state", "county", "commodity", "damagecause", "monthcode", "month", "acres", "loss")
#write.csv(RMA, file = "/dmine/code/git/paper-PHD1/data/RMA_csv/US_RMA_2001_2015.csv")
#save(RMA,file="/dmine/code/git/paper-PHD1/data/RMA_Rda/US_RMA_2001_2015.Rda")


RMA_PNW <- subset(RMA, state == "WA" | state == "OR" | state == "ID" )
#save(RMA_PNW,file="/dmine/code/git/paper-PHD1/data/RMA_Rda/PNW_RMA_2001_2015.Rda")
#write.csv(RMA_PNW, file = "/dmine/code/git/paper-PHD1/data/RMA_csv/RMA_PNW_2001_2015.csv")

RMA_PNW_1989 <- subset(RMA_1989, state == "WA" | state == "OR" | state == "ID" )
#save(RMA_PNW_1989,file="/dmine/code/git/paper-PHD1/data/RMA_Rda/PNW_RMA_1989_2000.Rda")
write.csv(RMA_PNW_1989, file = "/dmine/code/git/paper-PHD1/data/RMA_csv/RMA_PNW_1989_2000.csv")


RMA_PNW$lossperacre <- RMA_PNW$loss / RMA_PNW$acres

RMA_PNW <- subset(RMA_PNW, month != "")

RMA_PNW$cropyear <- RMA_PNW$year

RMA_PNW$cropyear[RMA_PNW$monthcode >= 10] <- RMA_PNW$year + 1

RMA_PNW <- subset(RMA_PNW, year != 2016)



RMA_commodity_loss <- aggregate(RMA$loss, by = list(RMA$year, RMA$state, RMA$county, RMA$commodity), FUN = sum)
colnames(RMA_commodity_loss) <- c("year", "state", "county", "commodity", "loss")
RMA_commodity_loss <- subset(RMA_commodity_loss, commodity != "ADJUSTED GROSS REVENUE")

RMA_commodity_count <- aggregate(RMA$loss, by = list(RMA$year, RMA$state, RMA$county, RMA$commodity), FUN = length)
colnames(RMA_commodity_count) <- c("year", "state", "county", "commodity", "count")
RMA_commodity_count <- subset(RMA_commodity_count, commodity != "ADJUSTED GROSS REVENUE")


RMA_commodity_acres <- aggregate(RMA$acres, by = list(RMA$year, RMA$state, RMA$county, RMA$commodity), FUN = sum)
colnames(RMA_commodity_acres) <- c("year", "state", "county", "commodity", "acres")
RMA_commodity_acres <- subset(RMA_commodity_acres, commodity != "ADJUSTED GROSS REVENUE")


RMA_damage_loss <- aggregate(RMA$loss, by = list(RMA$year, RMA$state, RMA$county, RMA$commodity, RMA$damagecause), FUN = sum)
colnames(RMA_damage_loss) <- c("year", "state", "county", "commodity", "damagecause", "loss")


RMA_damage_count <- aggregate(RMA$loss, by = list(RMA$year, RMA$state, RMA$county, RMA$commodity, RMA$damagecause), FUN = length)
colnames(RMA_damage_count) <- c("year", "state", "county", "commodity", "damagecause", "count")


RMA_damage_acres <- aggregate(RMA$acres, by = list(RMA$year, RMA$state, RMA$county, RMA$commodity, RMA$damagecause), FUN = sum)
colnames(RMA_damage_acres) <- c("year", "state", "county", "commodity", "damaecause", "acres")


RMA_commodity_combined <- cbind(RMA_commodity_loss, RMA_commodity_count$count, RMA_commodity_acres$acres)
colnames(RMA_commodity_combined) <- c("year", "state", "county", "commodity", "loss", "count", "acres")

RMA_commodity_combined$lossperacre <- RMA_commodity_combined$loss / RMA_commodity_combined$acres
RMA_commodity_combined$lossperclaim <- RMA_commodity_combined$loss / RMA_commodity_combined$count
RMA_commodity_combined$acresperclaim <- RMA_commodity_combined$acres / RMA_commodity_combined$count


RMA_damage_combined <- cbind(RMA_damage_loss, RMA_damage_count$count, RMA_damage_acres$acres)
colnames(RMA_damage_combined) <- c("year", "state", "county", "commodity", "damagecause", "loss", "count", "acres")
RMA_damage_combined$lossperacre <- RMA_damage_combined$loss / RMA_damage_combined$acres
RMA_damage_combined$lossperclaim <- RMA_damage_combined$loss / RMA_damage_combined$count
RMA_damage_combined$acresperclaim <- RMA_damage_combined$acres / RMA_damage_combined$count


RMA_damage_combined$lossperacre[!is.finite(RMA_damage_combined$lossperacre)] <- 0
#write.csv(RMA_commodity_combined, file = paste("/nfs/soilsesfeedback-data/data/RMA/RMA_commodity_combined.csv", sep=""))
#write.csv(RMA_damage_combined, file = paste("/nfs/soilsesfeedback-data/data/RMA/RMA_damage_combined.csv", sep=""))

#subset for PNW - 2001 - 2015
RMA_damage_PNW <- subset(RMA_damage_combined, state == "WA" | state == "OR" | state == "ID")

RMA_damage_PNW <- subset(RMA_damage_PNW, county != "All Other Counties")
#write.csv(RMA_damage_PNW, file = "/dmine/data/USDA/crop_indemnity_PNW/RMA_damage_PNW_2001_2015.csv")

RMA_Idaho <- subset(RMA_damage_PNW, state == "ID")
RMA_Oregon <- subset(RMA_damage_PNW, state == "OR")
RMA_Washington <- subset(RMA_damage_PNW, state == "WA")

RMA_Idaho_IPNW <- subset(RMA_Idaho, county == "Idaho" | county == "Lewis" | county ==  "Nez Perce" | county ==  "Latah" | county ==  "Benewah")
RMA_Oregon_IPNW <- subset(RMA_Oregon, county == "Wasco" | county ==  "Sherman" | county ==  "Gilliam" | county ==  "Morrow" | county ==  "Umatilla" | county ==  "Union" | county ==  "Wallowa")
RMA_Washington_IPNW <- subset(RMA_Washington, county ==  "Douglas" | county ==  "Grant" | county ==  "Benton" | county ==  "Franklin" | county ==  "Walla Walla" | county ==  "Adams" | county ==  "Lincoln" | county ==  "Spokane" | county ==  "Whitman" | county ==  "Columbia" | county ==  "Garfield" | county ==  "Asotin")

RMA_damage_IPNW <- rbind(RMA_Idaho_IPNW,RMA_Oregon_IPNW,RMA_Washington_IPNW)

#-load summarized data
Southern_ID_sumloss <- read.csv(text=getURL
("https://raw.githubusercontent.com/erichseamon/dmine_anova/master/PNW_summary_all.csv"), 
header = TRUE)

Southern_ID_sumloss_all_sum <- RMA_damage_PNW

Southern_ID_sumloss_all_sum  <- aggregate(loss ~ year + 
                                            damagecause + state + county + commodity,  
                                         Southern_ID_sumloss, sum)

Southern_ID_count_all_count  <- aggregate(count ~ year + 
                                            damagecause + county + commodity,  
                                          Southern_ID_sumloss, sum)
Southern_ID_sumloss_all_sum <- 
  Southern_ID_sumloss_all_sum[Southern_ID_sumloss_all_sum$loss >= 1, ]

#---subsetting for palouse/IPNW for sum, count, and lossperacre

#- SUM: sumloss_palouse_sum
#- COUNT: count_palouse_count
#- LOSSPERACRE: palouse_lossperacres

 palousecounties    <-  c("Idaho", "Lewis", "Nez Perce", "Clearwater", "Latah", "Benewah", "Kootenai","Douglas", "Grant", "Benton", "Franklin", "Walla Walla", "Adams", "Lincoln", "Spokane", "Whitman", "Columbia", "Garfield", "Asotin","Wasco", "Sherman", "Gilliam", "Morrow", "Umatilla", "Union", "Wallowa")

sumloss_palouse <- Southern_ID_sumloss[Southern_ID_sumloss$county %in% palousecounties, ]
sumloss_nonpalouse <- Southern_ID_sumloss[!Southern_ID_sumloss$county %in% palousecounties, ]
                      

sumloss_palouse_sum  <- aggregate(loss ~ year + 
                                            damagecause + county + state + commodity,  
                                          sumloss_palouse, sum)
count_palouse_count  <- aggregate(count ~ year + 
                                            damagecause + county + state + commodity,  
                                          sumloss_palouse, sum)

palouse_lossperacres  <- aggregate(lossperacres ~ year + 
                                            damagecause + county + state + commodity,  
                                          sumloss_palouse, mean)
sumloss_palouse_sum <- 
  sumloss_palouse_sum[sumloss_palouse_sum$loss >= 1, ]


#creating transforms for sum loss - cube root, square root, inverse, log10 and natural log, and then scale and center

Math.cbrt <- function(x) {
  sign(x) * abs(x)^(1/3)
}

#-calculate cube loss
sumloss_palouse_sum$cube_loss <- Math.cbrt(sumloss_palouse_sum$loss)

#-remove zeros
sumloss_palouse_sum <- subset(sumloss_palouse_sum, loss > 0)

#-use a log transform 

sumloss_palouse_sum$log10_loss <- log10(sumloss_palouse_sum$loss)
sumloss_palouse_sum$log_loss <- log(sumloss_palouse_sum$loss)

#-inverse transform
sumloss_palouse_sum$inverse_loss <- 1/sumloss_palouse_sum$loss

#sq root transformation
sumloss_palouse_sum$sqroot_loss <- sqrt(sumloss_palouse_sum$loss)

#--scale and center
sumloss_palouse_sum$scaled_cube_loss <- scale(sumloss_palouse_sum$cube_loss)
sumloss_palouse_sum$scaled_inverse_loss <- 
  scale(sumloss_palouse_sum$inverse_loss, center = TRUE, scale = TRUE)

#----


#-Loading all WHEAT claims for the palouse from 1989-2015
palouse_sumloss <- read.csv(text=getURL
("https://raw.githubusercontent.com/erichseamon/dmine_anova/master/Palouse_summary_sumloss.csv"), 
header = TRUE)
palouse_counts <- read.csv(text=getURL
("https://raw.githubusercontent.com/erichseamon/dmine_anova/master/Palouse_summary_counts.csv"), 
header = TRUE)


#use a cube transformation on loss for WHEAT claims
Math.cbrt <- function(x) {
  sign(x) * abs(x)^(1/3)
}

Southern_ID_sumloss_all_sum$cube_loss <- Math.cbrt(Southern_ID_sumloss_all_sum$loss)
Southern_ID_count_all_count$cube_counts <- Math.cbrt(Southern_ID_count_all_count$count)


#--aggregate palouse but not by county
palouse_sumloss_aggregate_county <- aggregate(loss ~ damagecause + year + commodity + county, 
                                       palouse_sumloss, mean)

#--aggregate palouse
palouse_sumloss_aggregate <- aggregate(loss ~ damagecause + year + commodity + county, 
                                       palouse_sumloss, sum)

#-calculate cube loss
palouse_sumloss_aggregate$cube_loss <- Math.cbrt(palouse_sumloss_aggregate$loss)

#-remove zeros
palouse_sumloss_aggregate <- subset(palouse_sumloss_aggregate, loss > 0)

#-use a log transform 

palouse_sumloss_aggregate$log10_loss <- log10(palouse_sumloss_aggregate$loss)
palouse_sumloss_aggregate$log_loss <- log(palouse_sumloss_aggregate$loss)

#-inverse transform
palouse_sumloss_aggregate$inverse_loss <- 1/palouse_sumloss_aggregate$loss

#sq root transformation
palouse_sumloss_aggregate$sqroot_loss <- sqrt(palouse_sumloss_aggregate$loss)

#--scale and center
palouse_sumloss_aggregate$scaled_cube_loss <- scale(palouse_sumloss_aggregate$cube_loss)
palouse_sumloss_aggregate$scaled_inverse_loss <- 
  scale(palouse_sumloss_aggregate$inverse_loss, center = TRUE, scale = TRUE)

#----


Step 2. Original Insurance Loss Dataset - Pacific Northwest

In Step 2, we document the following table, which is an example of the original data that was acquired from the USDA Risk Management Agency (RMA). Each record represents an individual insurance claim.

knitr::kable(RMA_PNW[1:10,], format="markdown") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
year state county commodity damagecause monthcode month acres loss lossperacre cropyear
11338 2001 ID Ada All Other Crops Drought 9 SEP 17.000 153.00 9.000000 2001
11339 2001 ID Ada All Other Crops Heat 8 AUG 105.200 5249.00 49.895437 2001
11340 2001 ID Ada All Other Crops Freeze 4 APR 125.000 4500.00 36.000000 2001
11341 2001 ID Ada All Other Crops Wind/Excess Wind 5 MAY 50.000 1800.00 36.000000 2001
11342 2001 ID Ada All Other Crops Wind/Excess Wind 4 APR 92.500 3330.00 36.000000 2001
11343 2001 ID Bannock WHEAT Drought 8 AUG 133.000 1212.00 9.112782 2001
11344 2001 ID Bannock WHEAT Drought 9 SEP 777.520 24807.00 31.905289 2001
11345 2001 ID Bannock WHEAT Drought 7 JUL 3529.754 54726.46 15.504327 2001
11346 2001 ID Bannock WHEAT Heat 7 JUL 19.796 2371.60 119.801980 2001
11347 2001 ID Bannock WHEAT Heat 8 AUG 25.000 904.00 36.160000 2001


Step 3. Aggregated Insurance Loss Dataset - Pacific Northwest

In Step 3, the following table is an example of an aggregated dataset derived from the original insurance loss csv files. Here we have summarized claims by year, county, commodity, and damage cause. Each unique combination is summarized, which includes the total summarized loss, the number of claims, the total summarized acreage, loss per acre, loss per claim, and acres per claim. This dataset was the basis for our data examination.

knitr::kable(RMA_damage_PNW[1:10,], format="markdown") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
year state county commodity damagecause loss count acres lossperacre lossperclaim acresperclaim
11 2014 ID Ada All Other Crops Area Plan Crops Only 1398 1 0 0.0000000 1398.0 0
12 2015 ID Ada All Other Crops Area Plan Crops Only 11810 1 0 0.0000000 11810.0 0
211 2008 OR Baker All Other Crops Area Plan Crops Only 15292 2 5482 2.7894929 7646.0 2741
212 2010 OR Baker All Other Crops Area Plan Crops Only 1819 2 2282 0.7971078 909.5 1141
215 2009 ID Bannock All Other Crops Area Plan Crops Only 2284 1 600 3.8066667 2284.0 600
245 2008 ID Bear Lake All Other Crops Area Plan Crops Only 8102 1 2468 3.2828201 8102.0 2468
246 2012 ID Bear Lake All Other Crops Area Plan Crops Only 7459 1 2468 3.0222853 7459.0 2468
291 2010 ID Bingham All Other Crops Area Plan Crops Only 4197 1 192 21.8593750 4197.0 192
292 2011 ID Bingham All Other Crops Area Plan Crops Only 7769 1 240 32.3708333 7769.0 240
293 2012 ID Bingham All Other Crops Area Plan Crops Only 6786 1 168 40.3928571 6786.0 168



Step 4: inland Pacific Northwest (iPNW) Study Area

In Step 4, we document the 24 county iPNW study area.

library(data.table)
library(maptools)
library(classInt)
library(leaflet)
library(dplyr)
library(raster)
library(htmltools)


states <- readShapePoly('/dmine/data/states/states.shp', 
                                  proj4string=CRS
                                  ("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
        projection = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
    

#setwd("/dmine/data/counties/")

counties <- readShapePoly('/dmine/data/counties/threestate_palouse.shp',
                     proj4string=CRS
                     ("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
projection = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")

#pal <- colorNumeric(palette = c("white", "orange", "darkorange", "red", "darkred"),
#                    domain = counties$NAME)
exte <- as.vector(extent(counties))

#label <- paste(sep = "<br/>", counties$NAME, round(counties$NAME, 0))
#markers <- data.frame(label)
#labs <- as.list(counties$NAME)

counties <- counties[counties$NAME != "Kootenai",]
counties <- counties[counties$NAME != "Clearwater",]

      lat_long <- coordinates(counties)

 counties_a <- data.frame(counties)
      labs <- lapply(seq(nrow(counties_a)), function(i) {
        paste0(as.character(counties_a[i, "NAME"]))
      })

#leaflet(data = counties, options = leafletOptions(minZoom = 6, maxZoom = 6)) %>%  addProviderTiles("Stamen.TonerLite") %>% fitBounds(exte[1], exte[3], exte[2], exte[4]) %>% addPolygons(color = "blue",  weight = 1) %>%

      
#--


mapz <- leaflet(data = counties, options = leafletOptions(zoomControl = FALSE, zoomSnap = .4, zoomDelta = .4)) %>%  
  addProviderTiles(providers$Hydda.Base) %>%
   addProviderTiles(providers$Stamen.TonerLines) %>% 
  
 
  
  fitBounds(exte[1], exte[3], exte[2], exte[4]) %>% addPolygons(color = "lightyellow",  fillOpacity = .8, weight = 1) %>%  addPolygons(data = states, fillOpacity = 0, weight = 4, stroke = TRUE, smoothFactor = 0.5, color = "black") %>%  addPolygons(color = "gray",  fillOpacity = 0, weight = 2) %>%  

  
 # addMiniMap(
#    tiles = providers$Stamen.TonerLite,
#    position = 'topleft', 
#    width = 175, height = 175,
#    toggleDisplay = FALSE) %>%
  
  
addLabelOnlyMarkers(data = counties, lng = lat_long[,1], lat = lat_long[,2], label = lapply(labs, HTML), labelOptions = labelOptions(style = list("font-family" = "serif"), noHide = TRUE, direction = 'middle', textOnly = TRUE, textsize = "16px", color = "white",  offset = c(20, 0), markerOptions(riseOnHover = TRUE, col = "white")
)) 
        
addScaleBar(mapz, position = c("topright"), options = scaleBarOptions())
#markerOptions(map, riseOnHover = TRUE)
 # addLegend(colors = ~NAME, labels = ~NAME, values = ~NAME, opacity = 0.5, title = NULL,
  #          position = "bottomright")



Step 5: PCA: PNW insurance loss, by COUNTY with COMMODITY loadings: 2001 to 2015

In Step 5, we perform a PCA for the insurance loss for the entire PNW, by county, with commodity as the factor loadings. Data from 2001 is 2015 is used. We additionally have generated a scree plot that shows the proportion of variance explained by the individual components.

#PNW by commodity

library(reshape2) ; RMA_damage_PNW_transformed_commodity <- dcast(RMA_damage_PNW, county + state ~ commodity, value.var = c("loss"), sum)
RMA_damage_PNW_transformed_commodity$statecounty <- paste(RMA_damage_PNW_transformed_commodity$county, "_", RMA_damage_PNW_transformed_commodity$state, sep="")
rownames(RMA_damage_PNW_transformed_commodity) <- RMA_damage_PNW_transformed_commodity$statecounty
RMA_damage_PNW_exogenous_commodity <- RMA_damage_PNW_transformed_commodity[,3:41]

RMA_damage_PNW_exogenous_commodity <- RMA_damage_PNW_exogenous_commodity[-1]
RMA_damage_PNW_exogenous_commodity <- RMA_damage_PNW_exogenous_commodity[-20]
RMA_damage_PNW_exogenous_commodity <- RMA_damage_PNW_exogenous_commodity[-37]


fit <- princomp(RMA_damage_PNW_exogenous_commodity, cor=TRUE)
#autoplot(prcomp(RMA_damage_PNW_exogenous_commodity, scale = TRUE, center = TRUE), data = RMA_damage_PNW_transformed_commodity, colour = 'county', loadings = TRUE, loadings.label = TRUE, loadings.label.size  = 3, label = FALSE, label.size = 3, legend = FALSE) + geom_text(vjust=-1, label=rownames(RMA_damage_PNW_exogenous_commodity)) + theme(legend.position = "none")

autoplot(prcomp(RMA_damage_PNW_exogenous_commodity, center = TRUE, scale = TRUE), frame = TRUE, frame.type = 'norm', data = RMA_damage_PNW_transformed_commodity, colour = 'state', loadings = TRUE, loadings.label = TRUE, loadings.label.size  = 3, label = FALSE, label.size = 2, legend = FALSE)  + theme_bw()   + geom_text(vjust=-1, label=rownames(RMA_damage_PNW_exogenous_commodity)) + theme(legend.position = "none") + ggtitle("Principle Components: PNW insurance loss\nby county with commodity loadings ") + theme(plot.title = element_text(size =28, family = "Serif")) + theme(axis.title.y = element_text(family = "Serif", size=16), axis.title.x = element_text(family = "Serif", size = 16), axis.text.x = element_text(size=rel(1.5), hjust = 1, family = "Serif"), axis.text.y = element_text(size=rel(1.5), hjust = 1, family = "Serif"))

#pca_RMA_damage_IPNW_exogenous$rotation

#plot(pca_RMA_damage_IPNW_exogenous$rotation)

std_dev <- fit$sdev

pr_var <- std_dev^2

#pr_var[1:10]

prop_varex <- pr_var/sum(pr_var)

plot(prop_varex, xlab = "Principal Component",
             ylab = "Proportion of Variance Explained",
             type = "b", main = "Scree Plot")


Step 6: PCA: PNW insurance loss, by YEAR with DAMAGE CAUSE loadings: 2001 to 2015

In Step 6, we perform a PCA for the insurance loss for the entire PNW, for all commodities by year, with damage cause as the factor loadings. Data from 2001 is 2015 is used. We additionally have generated a scree plot that shows the proportion of variance explained by the individual components.

  library(reshape2) ; RMA_damage_PNW_transformed_year <- dcast(RMA_damage_PNW, year ~ damagecause, value.var = c("loss"), sum)
#RMA_damage_PNW_transformed_year$statecounty <- paste(RMA_damage_PNW_transformed_year$county, "_", RMA_damage_PNW_transformed_year$state, sep="")
rownames(RMA_damage_PNW_transformed_year) <- RMA_damage_PNW_transformed_year$year
RMA_damage_PNW_exogenous_year <- RMA_damage_PNW_transformed_year[,2:31]

#PNW using year as group - with damage cause as loadings
#fit <- princomp(RMA_damage_PNW_exogenous_year, cor=TRUE, scale = TRUE, center = TRUE)

RMA_damage_PNW_exogenous_year_factor <- cbind(RMA_damage_PNW_exogenous_year[,2:3], RMA_damage_PNW_exogenous_year[6], RMA_damage_PNW_exogenous_year[8], RMA_damage_PNW_exogenous_year[,11:14], RMA_damage_PNW_exogenous_year[16:17])

fit <- princomp(RMA_damage_PNW_exogenous_year_factor, cor=TRUE)

autoplot(prcomp(RMA_damage_PNW_exogenous_year_factor, center = TRUE, scale = TRUE),  data = RMA_damage_PNW_exogenous_year_factor, colour = rownames(RMA_damage_PNW_exogenous_year_factor), loadings = TRUE, loadings.label = TRUE, loadings.label.size  = 3, label = FALSE, label.size = 2, legend = FALSE)  + theme_bw()   + geom_text(vjust=-1, label=rownames(RMA_damage_PNW_exogenous_year_factor)) + theme(legend.position = "none") + ggtitle("Principle Components: PNW insurance loss\nby year with damage cause loadings ") + theme(plot.title = element_text(size =28, family = "Serif")) + theme(axis.title.y = element_text(family = "Serif", size=16), axis.title.x = element_text(family = "Serif", size = 16), axis.text.x = element_text(size=rel(1.5), hjust = 1, family = "Serif"), axis.text.y = element_text(size=rel(1.5), hjust = 1, family = "Serif"))

#pca_RMA_damage_IPNW_exogenous$rotation

#plot(pca_RMA_damage_IPNW_exogenous$rotation)

std_dev <- fit$sdev

pr_var <- std_dev^2

#pr_var[1:10]

prop_varex <- pr_var/sum(pr_var)

plot(prop_varex, xlab = "Principal Component",
             ylab = "Proportion of Variance Explained",
             type = "b", main = "Scree Plot")


Step 7: PCA: PNW insurance loss, by MONTH with DAMAGE CAUSE loadings: 2001 to 2015

In Step 7, we perform a PCA for the insurance loss for the entire PNW, for all commodities by month, with damage cause as the factor loadings. Data from 2001 is 2015 is used. We additionally have generated a scree plot that shows the proportion of variance explained by the individual components.

#PNW damage cause as loadings using months as groups

#remove negatives related to lossperacre calc
RMA_PNW_noneg <- subset(RMA_PNW, lossperacre != Inf)
RMA_PNW_noneg2 <- subset(RMA_PNW_noneg, lossperacre >=0 )
RMA_PNW_noneg3 <- subset(RMA_PNW_noneg2, month !="" )

  library(reshape2) ; RMA_damage_PNW_transformed_month <- dcast(RMA_PNW_noneg3, month ~ damagecause, value.var = c("loss"), mean)

is.nan.data.frame <- function(x)
do.call(cbind, lapply(x, is.nan))

RMA_damage_PNW_transformed_month[is.nan(RMA_damage_PNW_transformed_month)] <- 0

#RMA_damage_PNW_transformed_year$statecounty <- paste(RMA_damage_PNW_transformed_year$county, "_", RMA_damage_PNW_transformed_year$state, sep="")
rownames(RMA_damage_PNW_transformed_month) <- RMA_damage_PNW_transformed_month$month
RMA_damage_PNW_exogenous_month <- RMA_damage_PNW_transformed_month[,2:31]

RMA_damage_PNW_exogenous_month <- cbind(RMA_damage_PNW_exogenous_month[,2:3], RMA_damage_PNW_exogenous_month[6], RMA_damage_PNW_exogenous_month[8], RMA_damage_PNW_exogenous_month[,11:14], RMA_damage_PNW_exogenous_month[16:17])


#autoplot(prcomp(RMA_damage_PNW_exogenous_year, center = TRUE, scale = TRUE), data = RMA_damage_PNW_transformed_year, colour = 'year', loadings = TRUE, loadings.label = TRUE, loadings.label.size  = 3, label = FALSE, label.size = 3, legend = FALSE) + geom_text(vjust=-1, label=rownames(RMA_damage_PNW_exogenous_year)) + theme(legend.position = "none")

#PNW using month as group - with damage cause as loadings
fit <- princomp(RMA_damage_PNW_exogenous_month, cor=TRUE)
autoplot(prcomp(RMA_damage_PNW_exogenous_month, center = TRUE, scale = TRUE),  data = RMA_damage_PNW_transformed_month, colour = 'month', loadings = TRUE, loadings.label = TRUE, loadings.label.size  = 3, label = FALSE, label.size = 2, legend = FALSE)  + theme_bw()   + geom_text(vjust=-1, label=rownames(RMA_damage_PNW_exogenous_month)) + theme(legend.position = "none") + ggtitle("Principle Components: PNW insurance loss\nby month with damage cause loadings ") + theme(plot.title = element_text(size =28, family = "Serif")) + theme(axis.title.y = element_text(family = "Serif", size=16), axis.title.x = element_text(family = "Serif", size = 16), axis.text.x = element_text(size=rel(1.5), hjust = 1, family = "Serif"), axis.text.y = element_text(size=rel(1.5), hjust = 1, family = "Serif"))

#pca_RMA_damage_IPNW_exogenous$rotation

#plot(pca_RMA_damage_IPNW_exogenous$rotation)

std_dev <- fit$sdev

pr_var <- std_dev^2

#pr_var[1:10]

prop_varex <- pr_var/sum(pr_var)

plot(prop_varex, xlab = "Principal Component",
             ylab = "Proportion of Variance Explained",
             type = "b", main = "Scree Plot")


Step 8: PCA: PNW WHEAT insurance loss, by YEAR with DAMAGE CAUSE loadings: 2001 to 2015

In Step 8, we perform a PCA for the insurance loss for the entire PNW, for wheat by year, with damage cause as the factor loadings. Data from 2001 is 2015 is used. We additionally have generated a scree plot that shows the proportion of variance explained by the individual components.

#PNW damage cause as loadings using years as groups, JUST FOR WHEAT

RMA_damage_PNW_wheat_year <- subset(RMA_damage_PNW, commodity == "WHEAT")
  library(reshape2) ; RMA_damage_PNW_transformed_wheat_year <- dcast(RMA_damage_PNW_wheat_year, year ~ damagecause, value.var = c("loss"), sum)
#RMA_damage_PNW_transformed_year$statecounty <- paste(RMA_damage_PNW_transformed_year$county, "_", RMA_damage_PNW_transformed_year$state, sep="")
rownames(RMA_damage_PNW_transformed_wheat_year) <- RMA_damage_PNW_transformed_wheat_year$year
RMA_damage_PNW_exogenous_wheat_year <- RMA_damage_PNW_transformed_wheat_year
par(mar=c(8, 6.1, 4, 2.1), family = 'serif', mgp=c(4, 1, 0), las=0)

#PNW using year as group - with damage cause as loadings: JUST FOR WHEAT

RMA_damage_PNW_exogenous_wheat_year_factor <- cbind(RMA_damage_PNW_exogenous_wheat_year[,3:4], RMA_damage_PNW_exogenous_wheat_year[7], RMA_damage_PNW_exogenous_wheat_year[9], RMA_damage_PNW_exogenous_wheat_year[,12:15], RMA_damage_PNW_exogenous_wheat_year[17:18])
fit <- princomp(RMA_damage_PNW_exogenous_wheat_year_factor, cor=TRUE)

autoplot(prcomp(RMA_damage_PNW_exogenous_wheat_year_factor, center = TRUE, scale = TRUE),  data = RMA_damage_PNW_exogenous_wheat_year_factor, colour = rownames(RMA_damage_PNW_exogenous_wheat_year_factor), loadings = TRUE, loadings.label = TRUE, loadings.label.size  = 3, label = FALSE, label.size = 2, legend = FALSE)  + theme_bw()   + geom_text(vjust=-1, label=rownames(RMA_damage_PNW_exogenous_wheat_year_factor)) + theme(legend.position = "none") + ggtitle("Principle Components: PNW WHEAT insurance loss\nby year with damage cause loadings ") + theme(plot.title = element_text(size =28, family = "Serif")) + theme(axis.title.y = element_text(family = "Serif", size=16), axis.title.x = element_text(family = "Serif", size = 16), axis.text.x = element_text(size=rel(1.5), hjust = 1, family = "Serif"), axis.text.y = element_text(size=rel(1.5), hjust = 1, family = "Serif"))

#pca_RMA_damage_IPNW_exogenous$rotation

#plot(pca_RMA_damage_IPNW_exogenous$rotation)

std_dev <- fit$sdev

pr_var <- std_dev^2

#pr_var[1:10]

prop_varex <- pr_var/sum(pr_var)

plot(prop_varex, xlab = "Principal Component",
             ylab = "Proportion of Variance Explained",
             type = "b", main = "Scree Plot")


Step 9: PCA: PNW WHEAT insurance loss, by COUNTY with DAMAGE CAUSE loadings: 2001 to 2015

In Step 9, we perform a PCA for the insurance loss for the entire PNW, for wheat by county, with damage cause as the factor loadings. Data from 2001 is 2015 is used. We additionally have generated a scree plot that shows the proportion of variance explained by the individual components.

#PNW wheat

RMA_damage_PNW_wheat <- subset(RMA_damage_PNW, commodity == "WHEAT")

library(reshape2) ; RMA_damage_PNW_transformed_wheat <- dcast(RMA_damage_PNW_wheat, county + state ~ damagecause, value.var = c("loss"), sum)
RMA_damage_PNW_transformed_wheat$statecounty <- paste(RMA_damage_PNW_transformed_wheat$county, "_", RMA_damage_PNW_transformed_wheat$state, sep="")
rownames(RMA_damage_PNW_transformed_wheat) <- RMA_damage_PNW_transformed_wheat$statecounty
#RMA_damage_PNW_exogenous_wheat <- RMA_damage_PNW_transformed_wheat[,3:27]
RMA_damage_PNW_exogenous_wheat <- RMA_damage_PNW_transformed_wheat

#autoplot(prcomp(RMA_damage_PNW_exogenous_month), data = RMA_damage_PNW_transformed_month, colour = 'month', loadings = TRUE, loadings.label = TRUE, loadings.label.size  = 3, label = FALSE, label.size = 3, legend = FALSE) + geom_text(vjust=-1, label=rownames(RMA_damage_PNW_exogenous_month)) + theme(legend.position = "none")

#PNW wheat

RMA_damage_PNW_exogenous_wheat <- cbind(RMA_damage_PNW_exogenous_wheat[,4:5], RMA_damage_PNW_exogenous_wheat[8], RMA_damage_PNW_exogenous_wheat[10], RMA_damage_PNW_exogenous_wheat[,13:16], RMA_damage_PNW_exogenous_wheat[18:19])

fit <- princomp(RMA_damage_PNW_exogenous_wheat, cor=TRUE)
#autoplot(prcomp(RMA_damage_PNW_exogenous_wheat, scale = TRUE, center = TRUE), data = RMA_damage_PNW_transformed_wheat, colour = 'county', loadings = TRUE, loadings.label = TRUE, loadings.label.size  = 3, label = FALSE, label.size = 3, legend = FALSE) + geom_text(vjust=-1, label=rownames(RMA_damage_PNW_exogenous_wheat)) + theme(legend.position = "none")

autoplot(prcomp(RMA_damage_PNW_exogenous_wheat, center = TRUE, scale = TRUE), frame = TRUE, frame.type = 'norm', data = RMA_damage_PNW_transformed_wheat, colour = 'state', loadings = TRUE, loadings.label = TRUE, loadings.label.size  = 3, label = FALSE, label.size = 2, legend = FALSE)  + theme_bw()   + geom_text(vjust=-1, label=rownames(RMA_damage_PNW_exogenous_wheat)) + theme(legend.position = "none") + ggtitle("Principle Components: PNW wheat insurance loss\nby county with damage cause loadings ") + theme(plot.title = element_text(size =28, family = "Serif")) + theme(axis.title.y = element_text(family = "Serif", size=16), axis.title.x = element_text(family = "Serif", size = 16), axis.text.x = element_text(size=rel(1.5), hjust = 1, family = "Serif"), axis.text.y = element_text(size=rel(1.5), hjust = 1, family = "Serif"))

#pca_RMA_damage_IPNW_exogenous$rotation

#plot(pca_RMA_damage_IPNW_exogenous$rotation)

std_dev <- fit$sdev

pr_var <- std_dev^2

#pr_var[1:10]

prop_varex <- pr_var/sum(pr_var)

plot(prop_varex, xlab = "Principal Component",
             ylab = "Proportion of Variance Explained",
             type = "b", main = "Scree Plot")


Step 10: PCA: PNW APPLES insurance loss, by COUNTY with DAMAGE CAUSE loadings: 2001 to 2015

In Step 5, we perform a PCA for the insurance loss for the entire PNW, for apples by county, with damage cause as the factor loadings. Data from 2001 is 2015 is used. We additionally have generated a scree plot that shows the proportion of variance explained by the individual components.

#PNW apples

RMA_damage_PNW_apples <- subset(RMA_damage_PNW, commodity == "APPLES")

  library(reshape2) ; RMA_damage_PNW_transformed_apples <- dcast(RMA_damage_PNW_apples, county + state ~ damagecause, value.var = c("loss"), sum)
RMA_damage_PNW_transformed_apples$statecounty <- paste(RMA_damage_PNW_transformed_apples$county, "_", RMA_damage_PNW_transformed_apples$state, sep="")
rownames(RMA_damage_PNW_transformed_apples) <- RMA_damage_PNW_transformed_apples$statecounty
RMA_damage_PNW_exogenous_apples <- RMA_damage_PNW_transformed_apples[,3:24]

RMA_damage_PNW_exogenous_apples <- cbind(RMA_damage_PNW_exogenous_apples[,1:4], RMA_damage_PNW_exogenous_apples[,7:10], RMA_damage_PNW_exogenous_apples[,13:14])

fit <- princomp(RMA_damage_PNW_exogenous_apples, cor=TRUE)
#PNW apples

#autoplot(prcomp(RMA_damage_PNW_exogenous_apples), data = RMA_damage_PNW_transformed_apples, colour = 'county', loadings = TRUE, loadings.label = TRUE, loadings.label.size  = 3, label = FALSE, label.size = 3, legend = FALSE) + geom_text(vjust=-1, label=rownames(RMA_damage_PNW_exogenous_apples)) + theme(legend.position = "none")

autoplot(prcomp(RMA_damage_PNW_exogenous_apples, center = TRUE, scale = TRUE), frame = TRUE, frame.type = 'norm', data = RMA_damage_PNW_transformed_apples, colour = 'state', loadings = TRUE, loadings.label = TRUE, loadings.label.size  = 3, label = FALSE, label.size = 2, legend = FALSE)  + theme_bw()   + geom_text(vjust=-1, label=rownames(RMA_damage_PNW_exogenous_apples)) + theme(legend.position = "none") + ggtitle("Principle Components: PNW apples insurance loss\nby county with damage cause loadings ") + theme(plot.title = element_text(size =28, family = "Serif")) + theme(axis.title.y = element_text(family = "Serif", size=16), axis.title.x = element_text(family = "Serif", size = 16), axis.text.x = element_text(size=rel(1.5), hjust = 1, family = "Serif"), axis.text.y = element_text(size=rel(1.5), hjust = 1, family = "Serif"))

#pca_RMA_damage_IPNW_exogenous$rotation

#plot(pca_RMA_damage_IPNW_exogenous$rotation)

std_dev <- fit$sdev

pr_var <- std_dev^2

#pr_var[1:10]

prop_varex <- pr_var/sum(pr_var)

plot(prop_varex, xlab = "Principal Component",
             ylab = "Proportion of Variance Explained",
             type = "b", main = "Scree Plot")


Step 11: PCA: PNW BARLEY insurance loss, by COUNTY with DAMAGE CAUSE loadings: 2001 to 2015

In Step 11, we perform a PCA for the insurance loss for the entire PNW, for barley by county, with damage cause as the factor loadings. Data from 2001 is 2015 is used. We additionally have generated a scree plot that shows the proportion of variance explained by the individual components.

#PNW barley

RMA_damage_PNW_barley <- subset(RMA_damage_PNW, commodity == "BARLEY")

  library(reshape2) ; RMA_damage_PNW_transformed_barley <- dcast(RMA_damage_PNW_barley, county + state ~ damagecause, value.var = c("loss"), sum)
RMA_damage_PNW_transformed_barley$statecounty <- paste(RMA_damage_PNW_transformed_barley$county, "_", RMA_damage_PNW_transformed_barley$state, sep="")
rownames(RMA_damage_PNW_transformed_barley) <- RMA_damage_PNW_transformed_barley$statecounty
#RMA_damage_PNW_exogenous_barley <- RMA_damage_PNW_transformed_barley[,3:23]
RMA_damage_PNW_exogenous_barley <- RMA_damage_PNW_transformed_barley

RMA_damage_PNW_exogenous_barley <- cbind(RMA_damage_PNW_exogenous_barley[,3:4], RMA_damage_PNW_exogenous_barley[7], RMA_damage_PNW_exogenous_barley[9], RMA_damage_PNW_exogenous_barley[,12:15], RMA_damage_PNW_exogenous_barley[17:18])

fit <- princomp(RMA_damage_PNW_exogenous_barley, cor=TRUE)

#PNW barley

#autoplot(prcomp(RMA_damage_PNW_exogenous_barley), data = RMA_damage_PNW_transformed_barley, colour = 'county', loadings = TRUE, loadings.label = TRUE, loadings.label.size  = 3, label = FALSE, label.size = 3, legend = FALSE) + geom_text(vjust=-1, label=rownames(RMA_damage_PNW_exogenous_barley)) + theme(legend.position = "none")

autoplot(prcomp(RMA_damage_PNW_exogenous_barley, center = TRUE, scale = TRUE), frame = TRUE, frame.type = 'norm', data = RMA_damage_PNW_transformed_barley, colour = 'state', loadings = TRUE, loadings.label = TRUE, loadings.label.size  = 3, label = FALSE, label.size = 2, legend = FALSE)  + theme_bw()   + geom_text(vjust=-1, label=rownames(RMA_damage_PNW_exogenous_barley)) + theme(legend.position = "none") + ggtitle("Principle Components: PNW barley insurance loss\nby county with damage cause loadings ") + theme(plot.title = element_text(size =28, family = "Serif")) + theme(axis.title.y = element_text(family = "Serif", size=16), axis.title.x = element_text(family = "Serif", size = 16), axis.text.x = element_text(size=rel(1.5), hjust = 1, family = "Serif"), axis.text.y = element_text(size=rel(1.5), hjust = 1, family = "Serif"))

#pca_RMA_damage_IPNW_exogenous$rotation

#plot(pca_RMA_damage_IPNW_exogenous$rotation)

std_dev <- fit$sdev

pr_var <- std_dev^2

#pr_var[1:10]

prop_varex <- pr_var/sum(pr_var)

plot(prop_varex, xlab = "Principal Component",
             ylab = "Proportion of Variance Explained",
             type = "b", main = "Scree Plot")


Step 12: PCA: IPNW insurance loss by COUNTY with DAMAGE CAUSE loadings: 2001 to 2015

In Step 5, we perform a PCA for the insurance loss for the inland PNW, for all commodities by county, with damage cause as the factor loadings. Data from 2001 is 2015 is used. We additionally have generated a scree plot that shows the proportion of variance explained by the individual components.

#all PNW by damage

  library(reshape2) ; RMA_damage_PNW_transformed <- dcast(RMA_damage_PNW, county + state ~ damagecause, value.var = c("loss"), sum)
RMA_damage_PNW_transformed$statecounty <- paste(RMA_damage_PNW_transformed$county, "_", RMA_damage_PNW_transformed$state, sep="")
rownames(RMA_damage_PNW_transformed) <- RMA_damage_PNW_transformed$statecounty
RMA_damage_PNW_exogenous <- RMA_damage_PNW_transformed[,3:32]
par(mar=c(8, 6.1, 4, 2.1), family = 'serif', mgp=c(4, 1, 0), las=0)

fit <- stats::princomp(RMA_damage_PNW_exogenous, cor=TRUE)
autoplot(prcomp(RMA_damage_PNW_exogenous, center = TRUE, scale = TRUE), frame = TRUE, frame.type = 'norm', data = RMA_damage_PNW_transformed, colour = 'state', loadings = TRUE, loadings.label = TRUE, loadings.label.size  = 3, label = FALSE, label.size = 2, legend = FALSE)  + theme_bw()   + geom_text(vjust=-1, label=rownames(RMA_damage_PNW_exogenous)) + theme(legend.position = "none") + ggtitle("Principle Components: IPNW insurance loss\nby county with damage cause loadings: 2001 to 2015") + theme(plot.title = element_text(size =28, family = "Serif")) + theme(axis.title.y = element_text(family = "Serif", size=16), axis.title.x = element_text(family = "Serif", size = 16), axis.text.x = element_text(size=rel(1.5), hjust = 1, family = "Serif"), axis.text.y = element_text(size=rel(1.5), hjust = 1, family = "Serif"))

par(mar=c(8, 6.1, 4, 2.1), family = 'serif', mgp=c(4, 1, 0), las=0)
#pca_RMA_damage_IPNW_exogenous$rotation

#plot(pca_RMA_damage_IPNW_exogenous$rotation)

std_dev <- fit$sdev

pr_var <- std_dev^2

#pr_var[1:10]

prop_varex <- pr_var/sum(pr_var)

plot(prop_varex, xlab = "Principal Component",
             ylab = "Proportion of Variance Explained",
             type = "b", main = "Scree Plot")


Step 13: PCA: IPNW insurance loss, by COUNTY with DAMAGE CAUSE loadings: 2001 to 2015

In Step 5, we perform a PCA for the insurance loss for the inland PNW, for all commodities by county, with damage cause as the factor loadings. Data from 2001 is 2015 is used. We additionally have generated a scree plot that shows the proportion of variance explained by the individual components.

#just the IPNW

RMA_damage_IPNW$log10loss <- log10(RMA_damage_IPNW$loss)


RMA_damage_IPNW_wheat <- subset(RMA_damage_IPNW, commodity == "WHEAT")
 library(reshape2) ; RMA_damage_IPNW_transformed <- dcast(RMA_damage_IPNW_wheat, county + state ~ damagecause, value.var = c("loss"), sum)
RMA_damage_IPNW_transformed$statecounty <- paste(RMA_damage_IPNW_transformed$county, "_", RMA_damage_IPNW_transformed$state, sep="")
rownames(RMA_damage_IPNW_transformed) <- RMA_damage_IPNW_transformed$statecounty
RMA_damage_IPNW_exogenous <- RMA_damage_IPNW_transformed[,3:26]


#IPNW

#RMA_damage_IPNW_exogenous <- scale(RMA_damage_IPNW_exogenous, scale = TRUE, center = TRUE)
pca_RMA_damage_IPNW_exogenous <- prcomp(RMA_damage_IPNW_exogenous, scale = TRUE, center = TRUE)

RMA_damage_IPNW_exogenous_factor <- cbind(RMA_damage_IPNW_exogenous[2:3], RMA_damage_IPNW_exogenous[5:6], RMA_damage_IPNW_exogenous[8], RMA_damage_IPNW_exogenous[11], RMA_damage_IPNW_exogenous[13:14], RMA_damage_IPNW_exogenous[16:17])

fit <- princomp(RMA_damage_IPNW_exogenous_factor, cor=TRUE)

autoplot(prcomp(RMA_damage_IPNW_exogenous_factor, scale = TRUE, center = TRUE), frame = TRUE, frame.type = 'norm', data = RMA_damage_IPNW_transformed, colour = 'state', loadings = TRUE, loadings.label = TRUE, loadings.label.size  = 3, label = FALSE, label.size = 2, legend = FALSE)  + theme_bw()   + geom_text(vjust=-1, label=rownames(RMA_damage_IPNW_exogenous)) + theme(legend.position = "none") + ggtitle("Principle Components: IPNW insurance loss\nby county with damage cause loadings ") + theme(plot.title = element_text(size =28, family = "Serif")) + theme(axis.title.y = element_text(family = "Serif", size=16), axis.title.x = element_text(family = "Serif", size = 16), axis.text.x = element_text(size=rel(1.5), hjust = 1, family = "Serif"), axis.text.y = element_text(size=rel(1.5), hjust = 1, family = "Serif"))

#pca_RMA_damage_IPNW_exogenous$rotation

#plot(pca_RMA_damage_IPNW_exogenous$rotation)

std_dev <- fit$sdev

pr_var <- std_dev^2

#pr_var[1:10]

prop_varex <- pr_var/sum(pr_var)

plot(prop_varex, xlab = "Principal Component",
             ylab = "Proportion of Variance Explained",
             type = "b", main = "Scree Plot")


Step 14: PCA: IPNW insurance loss, by YEAR with DAMAGE CAUSE loadings: 2001 to 2015

In Step 5, we perform a PCA for the insurance loss for the inland PNW, for all commodities by year, with damage cause as the factor loadings. Data from 2001 is 2015 is used. We additionally have generated a scree plot that shows the proportion of variance explained by the individual components.

#IPNW damage cause as loadings using years as groups

 library(reshape2) ; RMA_damage_IPNW_transformed_year <- dcast(RMA_damage_IPNW, year ~ damagecause, value.var = c("loss"), sum)
#RMA_damage_PNW_transformed_year$statecounty <- paste(RMA_damage_PNW_transformed_year$county, "_", RMA_damage_PNW_transformed_year$state, sep="")
rownames(RMA_damage_IPNW_transformed_year) <- RMA_damage_IPNW_transformed_year$year
#RMA_damage_IPNW_exogenous_year <- RMA_damage_IPNW_transformed_year[,2:28]
RMA_damage_IPNW_exogenous_year <- RMA_damage_IPNW_transformed_year

#IPNW with years as the group

RMA_damage_IPNW_exogenous_year <- cbind(RMA_damage_IPNW_exogenous_year[3:4], RMA_damage_IPNW_exogenous_year[7], RMA_damage_IPNW_exogenous_year[9], RMA_damage_IPNW_exogenous_year[,12:15], RMA_damage_IPNW_exogenous_year[18:19])

fit <- princomp(RMA_damage_IPNW_exogenous_year, cor=TRUE)

autoplot(prcomp(RMA_damage_IPNW_exogenous_year, center = TRUE, scale = TRUE), data = RMA_damage_IPNW_transformed_year, loadings = TRUE, loadings.label = TRUE, loadings.label.size  = 3, label = FALSE, label.size = 2, legend = FALSE)  + theme_bw()  + geom_text(vjust=-1, label=rownames(RMA_damage_IPNW_exogenous_year)) + theme(legend.position = "none") + ggtitle("Principle Components: IPNW insurance loss\nby year with damage cause loadings") + theme(plot.title = element_text(size =28, family = "Serif")) + theme(axis.title.y = element_text(family = "Serif", size=16), axis.title.x = element_text(family = "Serif", size = 16), axis.text.x = element_text(size=rel(1.5), hjust = 1, family = "Serif"), axis.text.y = element_text(size=rel(1.5), hjust = 1, family = "Serif"))

#autoplot(prcomp(RMA_damage_IPNW_exogenous, scale = TRUE, center = TRUE), data = RMA_damage_IPNW_transformed, colour = 'county', loadings = TRUE, loadings.label = TRUE, loadings.label.size  = 3, label = FALSE, label.size = 3, legend = FALSE) + geom_text(vjust=-1, label=rownames(RMA_damage_IPNW_exogenous)) + theme(legend.position = "none")
#pca_RMA_damage_IPNW_exogenous$rotation

#plot(pca_RMA_damage_IPNW_exogenous$rotation)

std_dev <- fit$sdev

pr_var <- std_dev^2

#pr_var[1:10]

prop_varex <- pr_var/sum(pr_var)

plot(prop_varex, xlab = "Principal Component",
             ylab = "Proportion of Variance Explained",
             type = "b", main = "Scree Plot")


Step 15: PCA: IPNW insurance loss, by COUNTY with DAMAGE CAUSE loadings: 2015

In Step 5, we perform a PCA for the insurance loss for the inland PNW, for all commodities by year, with damage cause as the factor loadings. Data for just 2015 is used. We additionally have generated a scree plot that shows the proportion of variance explained by the individual components.

#just the IPNW 2015

RMA_damage_IPNW_2015 <- subset(RMA_damage_IPNW, year == 2009)

 library(reshape2) ; RMA_damage_IPNW_transformed_2015 <- dcast(RMA_damage_IPNW_2015, county + state ~ damagecause, value.var = c("loss"), sum)
RMA_damage_IPNW_transformed_2015$statecounty <- paste(RMA_damage_IPNW_transformed_2015$county, "_", RMA_damage_IPNW_transformed_2015$state, sep="")
rownames(RMA_damage_IPNW_transformed_2015) <- RMA_damage_IPNW_transformed_2015$statecounty
#RMA_damage_IPNW_exogenous_2015 <- RMA_damage_IPNW_transformed_2015[,3:21]
RMA_damage_IPNW_exogenous_2015 <- RMA_damage_IPNW_transformed_2015

#IPNW 2015
RMA_damage_IPNW_exogenous_2015 <- cbind(RMA_damage_IPNW_exogenous_2015[4:5], RMA_damage_IPNW_exogenous_2015[,8:13], RMA_damage_IPNW_exogenous_2015[,15:16])

fit <- princomp(RMA_damage_IPNW_exogenous_year, cor=TRUE)

#RMA_damage_IPNW_exogenous <- scale(RMA_damage_IPNW_exogenous, scale = TRUE, center = TRUE)
pca_RMA_damage_IPNW_exogenous_2015 <- prcomp(RMA_damage_IPNW_exogenous_2015, scale = TRUE, center = TRUE)
#autoplot(prcomp(RMA_damage_IPNW_exogenous_2015, scale = TRUE, center = TRUE), data = RMA_damage_IPNW_transformed_2015, colour = 'state', loadings = TRUE, loadings.label = TRUE, loadings.label.size  = 3, label = FALSE, label.size = 3, legend = FALSE) + geom_text(vjust=-1, label=rownames(RMA_damage_IPNW_exogenous_2015)) + theme(legend.position = "none")

autoplot(prcomp(RMA_damage_IPNW_exogenous_2015, center = TRUE, scale = TRUE), frame = TRUE, frame.type = 'norm', data = RMA_damage_IPNW_transformed_2015, colour = 'state', loadings = TRUE, loadings.label = TRUE, loadings.label.size  = 3, label = FALSE, label.size = 2, legend = FALSE)  + theme_bw()   + geom_text(vjust=-1, label=rownames(RMA_damage_IPNW_exogenous_2015)) + theme(legend.position = "none") + ggtitle("Principle Components: 2015 IPNW insurance loss\nby county with damage cause loadings") + theme(plot.title = element_text(size =28, family = "Serif")) + theme(axis.title.y = element_text(family = "Serif", size=16), axis.title.x = element_text(family = "Serif", size = 16), axis.text.x = element_text(size=rel(1.5), hjust = 1, family = "Serif"), axis.text.y = element_text(size=rel(1.5), hjust = 1, family = "Serif"))

#pca_RMA_damage_IPNW_exogenous$rotation

#plot(pca_RMA_damage_IPNW_exogenous$rotation)

std_dev <- fit$sdev

pr_var <- std_dev^2

#pr_var[1:10]

prop_varex <- pr_var/sum(pr_var)

plot(prop_varex, xlab = "Principal Component",
             ylab = "Proportion of Variance Explained",
             type = "b", main = "Scree Plot")


Step 16: PCA: IPNW WHEAT loss by YEAR with DAMAGE CAUSE loadings: Temporal Mapping of PC1 and PC2

In Step 16, we summarize our PCA findings for the inland PNW, for wheat by year, using damage cause as the factor loadings. We take PC1 and PC2 loadings and plot those for each year.

library(htmlTable)

#year factanal

#factanal
#fit <- factanal(x = RMA_damage_PNW_exogenous_wheat_year_factor, factors = 3, rotation = "varimax")
#load <- fit$loadings[,1:2] 
#plot(load,type="n") # set up plot 
#text(load,labels=names(RMA_damage_PNW_exogenous_wheat_year_factor),cex=.7) # add variable names

library(FactoMineR)
result <- FactoMineR::PCA(RMA_damage_PNW_exogenous_wheat_year_factor) # graphs generated automatically

#--create table of PC1 and PC2 loadings

PC_table <- cbind(result$var$coord[,1], result$var$coord[,2])
colnames(PC_table) <- c("PC1", "PC2")
PC_table <- as.data.frame(PC_table)
PC_table$PC1 <- round(PC_table$PC1, 3)
PC_table$PC2 <- round(PC_table$PC2, 3)

htmlTable(PC_table,
          cgroup = c("2001-2015"),
          n.cgroup = c(2),
          caption = "Damage cause variable loadings by year",
          align = "c",
          rnames = rownames(PC_table),
          css.cell = "padding-left: .5em; padding-right: .5em;")
Damage cause variable loadings by year
2001-2015
PC1 PC2
Cold Wet Weather -0.325 0.919
Cold Winter 0.834 0.462
Drought 0.845 0.111
Excess Moisture/Precip/Rain -0.174 0.961
Fire 0.588 0.123
Flood -0.431 0.832
Freeze 0.729 0.396
Frost 0.168 -0.169
Heat 0.865 0.008
Hot Wind 0.365 0.015
#--year end

PC_table_year <- cbind(result$ind$coord[,1], result$ind$coord[,2])
colnames(PC_table_year) <- c("PC1", "PC2")

par(mar=c(6, 6.1, 4, 2.1), family = 'serif', mgp=c(4, 1, 0), las=0)

barplot(t(PC_table_year), beside = T, las = 2, col=c("blue", "red"), cex.axis = 1.5, cex.names = 1.5)

#legend("top", c("PC1","PC2"), pch=15, 
#col=c("blue","red"), 
#bty="n")


Step 17: PCA: IPNW WHEAT loss by COUNTY with DAMAGE CAUSE loadings: Spatial Mapping of PC1 and PC2

In Step 17, we summarize our PCA findings for the inland PNW, for wheat by county, using damage cause as the factor loadings. We take PC1 and PC2 loadings and plot each as a map.

#factanal
#fit <- factanal(x = RMA_damage_IPNW_exogenous_factor, factors = 3, rotation = "varimax")
#load <- fit$loadings[,1:2] 
#plot(load,type="n") # set up plot 
#text(load,labels=names(RMA_damage_IPNW_exogenous_factor),cex=.7) # add variable names

library(FactoMineR)
result <- FactoMineR::PCA(RMA_damage_IPNW_exogenous_factor) # graphs generated automatically

#--create table of PC1 and PC2 loadings

PC_table <- cbind(result$var$coord[,1], result$var$coord[,2])
colnames(PC_table) <- c("PC1", "PC2")
PC_table <- as.data.frame(PC_table)
PC_table$PC1 <- round(PC_table$PC1, 3)
PC_table$PC2 <- round(PC_table$PC2, 3)

htmlTable(PC_table,
          cgroup = c("2001-2015"),
          n.cgroup = c(2),
          caption = "Damage cause variable loadings by county",
          align = "c",
          rnames = rownames(PC_table),
          css.cell = "padding-left: .5em; padding-right: .5em;")
Damage cause variable loadings by county
2001-2015
PC1 PC2
Cold Wet Weather 0.465 0.797
Cold Winter 0.58 -0.244
Decline in Price 0.579 0.56
Drought 0.86 -0.21
Excess Moisture/Precip/Rain -0.107 0.897
Fire 0.738 0.043
Freeze 0.913 -0.176
Frost 0.945 -0.083
Heat 0.929 0.079
Hot Wind 0.882 -0.163
#-map PC loadings

PC1_map <- as.data.frame(result$ind$coord[,1])
PC1_map$NAME <- rownames(PC1_map)
colnames(PC1_map) <- c("PC1", "NAME")
PC1_NAME <- cbind(PC1_map, t(as.data.frame(strsplit(PC1_map$NAME, split='_', fixed=TRUE))))
colnames(PC1_NAME) <- c("PC1", "Statecounty", "NAME", "STATE")

PC1_map <- merge(counties, PC1_NAME, by=c("NAME"), all = FALSE)

PC2_map <- as.data.frame(result$ind$coord[,2])
PC2_map$NAME <- rownames(PC2_map)
colnames(PC2_map) <- c("PC2", "NAME")
PC2_NAME <- cbind(PC2_map, t(as.data.frame(strsplit(PC2_map$NAME, split='_', fixed=TRUE))))
colnames(PC2_NAME) <- c("PC2", "Statecounty", "NAME", "STATE")

PC2_map <- merge(counties, PC2_NAME, by=c("NAME"), all = FALSE)


#PC1 map

#-----

PC1_map$PC1[is.na(PC1_map$PC1)] <- 0

 
## Make vector of colors for values smaller than 0 (20 colors)
rc1 <- colorRampPalette(colors = c("blue", "lightblue", "white"), space = "Lab")(4)

## Make vector of colors for values larger than 0 (180 colors)
rc2 <- colorRampPalette(colors = c("white", "red"), space = "Lab")(16)

## Combine the two color palettes
rampcols <- c(rc1, rc2)



# pal <- colorNumeric(colorRampPalette(c("blue", "lightblue", "white", "lightpink", "red"))(7), na.color = "#ffffff",domain = eval(parse(text=paste("PC1_map$", "PC1", sep=""))))
 pal <- colorNumeric(palette = rampcols, na.color = "#ffffff",domain = eval(parse(text=paste("PC1_map$", "PC1", sep=""))))

map <- leaflet(data = PC1_map, options = leafletOptions(zoomControl = FALSE, zoomSnap = 0.1)) %>%  
  addProviderTiles(providers$Hydda.Base) %>%
   addProviderTiles(providers$Stamen.TonerLines) %>%
  
  addMiniMap(
    tiles = providers$Stamen.TonerLite,
    position = 'topleft', 
    width = 100, height = 100,
    toggleDisplay = FALSE) %>%
  
  
   #addControl(title, position = "topleft", className="map-title") %>%
  
  
  fitBounds(exte[1], exte[3], exte[2], exte[4]) %>% addPolygons(fillColor = ~pal(PC1_map$PC1),  fillOpacity = .8, weight = 1, stroke = TRUE, smoothFactor = 0.5, color = "black") %>% addPolygons(data = states, fillOpacity = 0, weight = 4, stroke = TRUE, smoothFactor = 0.5, color = "black") %>% 
  
    addLegend(pal = pal, values = na.omit(~PC1), group = "circles", title = paste("PC1 Loadings", sep = ""), na.label = "", position = "topright") %>%

  
  #label = lapply(labs, HTML),

addLabelOnlyMarkers(data = PC1_map, lng = lat_long[,1], lat = lat_long[,2],  label = lapply(round(PC1_map$PC1, 2), HTML), labelOptions = labelOptions(style = list("font-family" = "serif"), noHide = TRUE, direction = 'middle', textOnly = TRUE, textsize = "20px", col = "white",  offset = c(20, 0), markerOptions(riseOnHover = TRUE)
)) 
        
addScaleBar(map, position = c("topright"), options = scaleBarOptions())
#PC2 map

 
PC2_map$PC2[is.na(PC2_map$PC2)] <- 0

 
## Make vector of colors for values smaller than 0 (20 colors)
rc1 <- colorRampPalette(colors = c("blue", "lightblue", "white"), space = "Lab")(4)

## Make vector of colors for values larger than 0 (180 colors)
rc2 <- colorRampPalette(colors = c("white", "red"), space = "Lab")(16)

## Combine the two color palettes
rampcols <- c(rc1, rc2)

#-----
 #pal <- colorNumeric(brewer.pal(11,"RdBu"), na.color = "#ffffff",domain = eval(parse(text=paste("PC2_map$", "PC2", sep=""))))
 pal <- colorNumeric(palette = rampcols, na.color = "#ffffff",domain = eval(parse(text=paste("PC2_map$", "PC2", sep=""))))

map <- leaflet(data = PC2_map, options = leafletOptions(zoomControl = FALSE, zoomSnap = 0.1)) %>%  
  addProviderTiles(providers$Hydda.Base) %>%
   addProviderTiles(providers$Stamen.TonerLines) %>%
  
  addMiniMap(
    tiles = providers$Stamen.TonerLite,
    position = 'topleft', 
    width = 100, height = 100,
    toggleDisplay = FALSE) %>%
  
  
   #addControl(title, position = "topleft", className="map-title") %>%
  
  
  fitBounds(exte[1], exte[3], exte[2], exte[4]) %>% addPolygons(fillColor = ~pal(PC2_map$PC2),  fillOpacity = .8, weight = 1, stroke = TRUE, smoothFactor = 0.5, color = "black") %>% addPolygons(data = states, fillOpacity = 0, weight = 4, stroke = TRUE, smoothFactor = 0.5, color = "black") %>% 
  
    addLegend(pal = pal, values = na.omit(~PC2), group = "circles", title = paste("PC2 Loadings", sep = ""), na.label = "", position = "topright") %>%

  
  

addLabelOnlyMarkers(data = PC2_map, lng = lat_long[,1], lat = lat_long[,2], label = lapply(round(PC2_map$PC2, 2), HTML), labelOptions = labelOptions(style = list("font-family" = "serif"), noHide = TRUE, direction = 'middle', textOnly = TRUE, textsize = "20px", col = "white",  offset = c(20, 0), markerOptions(riseOnHover = TRUE)
)) 
        
addScaleBar(map, position = c("topright"), options = scaleBarOptions())